set.seed(42)
library(ggplot2) # for autoplot() generic
library(dplyr)
library(sf)
library(spgwr) ## Maybe change to GWmodel for parallel proc.
library(tmap)
dat_sf <- st_read("./data/glacier_clim.shp")
## Reading layer `glacier_clim' from data source `/Users/u0784726/Dropbox/Data/devtools/glacier_mb/data/glacier_clim.shp' using driver `ESRI Shapefile'
## Simple feature collection with 95086 features and 25 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 67.47845 ymin: 27.4918 xmax: 103.8915 ymax: 45.35089
## Geodetic CRS: GCS_unknown
dat_sf <- dat_sf[sample(1:nrow(dat_sf), 5000), ]
dat_sp <- as_Spatial(dat_sf)
Model formula
f1<- mb_mwea ~ t2m_d + tp_d +
z_med + z_slope + z_aspct +
area_m2 + tau
Bandwidth selection
bw <- gwr.sel(f1, dat_sp, adapt = TRUE)
## Adaptive q: 0.381966 CV score: 592.3786
## Adaptive q: 0.618034 CV score: 609.4785
## Adaptive q: 0.236068 CV score: 572.0944
## Adaptive q: 0.145898 CV score: 549.3333
## Adaptive q: 0.09016994 CV score: 531.1615
## Adaptive q: 0.05572809 CV score: 518.017
## Adaptive q: 0.03444185 CV score: 505.9349
## Adaptive q: 0.02128624 CV score: 494.9289
## Adaptive q: 0.01315562 CV score: 490.5888
## Adaptive q: 0.008130619 CV score: 490.0614
## Adaptive q: 0.009033444 CV score: 490.16
## Adaptive q: 0.005024999 CV score: 515.8404
## Adaptive q: 0.006944377 CV score: 492.4375
## Adaptive q: 0.007677515 CV score: 491.2669
## Adaptive q: 0.008475467 CV score: 490.0292
## Adaptive q: 0.00843161 CV score: 490.0452
## Adaptive q: 0.008688595 CV score: 490.1721
## Adaptive q: 0.008556875 CV score: 490.0597
## Adaptive q: 0.008516157 CV score: 490.0364
## Adaptive q: 0.008475467 CV score: 490.0292
print(bw)
## [1] 0.008475467
Run model
mb_gwr <- gwr(f1, dat_sp, adapt = bw)
## Warning in proj4string(data): CRS object has comment, which is lost in output
summary(mb_gwr)
## Length Class Mode
## SDF 5000 SpatialPointsDataFrame S4
## lhat 1 -none- logical
## lm 11 -none- list
## results 0 -none- NULL
## bandwidth 5000 -none- numeric
## adapt 1 -none- numeric
## hatmatrix 1 -none- logical
## gweight 1 -none- character
## gTSS 1 -none- numeric
## this.call 4 -none- call
## fp.given 1 -none- logical
## timings 8 -none- numeric
tm_shape(mb_gwr$SDF) +
tm_symbols(col = "localR2", size = 0.25, alpha = 0.75, style = "quantile")
## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
tm_shape(mb_gwr$SDF) +
tm_symbols(col = "t2m_d", size = 0.25, alpha = 0.75, style = "quantile")
## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
## Variable(s) "t2m_d" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
tm_shape(mb_gwr$SDF) +
tm_symbols(col = "tp_d", size = 0.25, alpha = 0.75, style = "quantile")
## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
## Variable(s) "tp_d" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
tm_shape(mb_gwr$SDF) +
tm_symbols(col = "z_med", size = 0.25, alpha = 0.75, style = "quantile")
## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
## Variable(s) "z_med" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
tm_shape(mb_gwr$SDF) +
tm_symbols(col = "tau", size = 0.25, alpha = 0.75, style = "quantile")
## Warning in sp::proj4string(obj): CRS object has comment, which is lost in output
## Variable(s) "tau" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.